% In this practical we will work with a single subject's data from an
% emotional faces task (data courtesy of Susie Murphy) and perform an
% time-frequency analysis in sensor space. 
% This dataset can be downloaded from:  
%  
% www.fmrib.ox.ac.uk/~woolrich/faces_subject1_data.tar.gz
% 
% Note that this contains the spm file:
% spm8_meg1.mat
% that is an SPM MEEG object that has continuous data that has already been
% SSS Maxfiltered and downsampled to 250 Hz. 
% 
% and
% espm8_meg1.mat
%
% which is an SPM MEEG object that has the same data epoched into the
% different task conditions.
%

%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

global OSLDIR;
    
% set this to where you have downloaded OSL and the practical data:
practical_dir='/home/mwoolrich/Desktop'; 
osldir=[practical_dir '/osl1.3.1'];    

%practical_dir='/Users/woolrich';
%osldir=[practical_dir '/homedir/matlab/osl1.3.1'];    

addpath(osldir);
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%
%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS

workingdir=[practical_dir '/faces_subject1_data']; % directory where the data is
workingdir=[practical_dir '/homedir/matlab/osl_testdata_dir/faces_subject1_data']; % directory where the data is

cmd = ['mkdir ' workingdir]; if ~exist(workingdir, 'dir'), unix(cmd); end% make dir to put the results in

clear spm_files spm_files_epoched;
% set up a list of SPM MEEG object file names (we only have one here)
spm_files{1}=[workingdir '/spm8_meg1.mat'];
spm_files_epoched{1}=[workingdir '/espm8_meg1.mat'];

%%%%%%%%%%%%%%%%%%%
%% SETUP SENSOR SPACE MULTI-BAND TIME-FREQ OAT
 
oat=[];
oat.source_recon.D_continuous=spm_files;
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
oat.source_recon.D_epoched=spm_files_epoched; % this is passed in so that the bad trials and bad channels can be read out
oat.source_recon.freq_range=[1 40]; % frequency range in Hz
oat.source_recon.time_range=[-0.4 0.5];
oat.source_recon.method='none';
oat.source_recon.dirname=[oat.source_recon.D_epoched{1} '_hilbert_multiband'];

% TF settings:
oat.first_level.tf_method='morlet'; % can be morlet or hilbert
oat.first_level.tf_freq_range=[5 40]; % frequency range in Hz
oat.first_level.time_range=[-0.2 0.3]; % need to make this time range smallet than oat.source_recon.time_range to remove edge effects
oat.first_level.tf_num_freqs=14; % we are keeping this unusally low in the practical for the sake of speed
%oat.first_level.tf_hilbert_freq_res=8;

% NOTE that you can also set HILBERT freq ranges explicitly, e.g.:
% oat.first_level.tf_hilbert_freq_ranges=[[4 8];[8 12];[12 16];[16 20];[20 24];[24 30]]; % frequency range in Hz

oat.first_level.post_tf_downsample_factor=4; % does downsampling after the TF decomposition

oat.first_level.bc=[1 1 0]; % specifies whether or not baseline correction is done for the different contrasts

% Xsummary is a parsimonious description of the design matrix.
% It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to croat_settingse the (num_regressors x num_trials) design matrix:
Xsummary={};
Xsummary{1}=[1 0 0 0];Xsummary{2}=[0 1 0 0];Xsummary{3}=[0 0 1 0];Xsummary{4}=[0 0 0 1];
oat.first_level.design_matrix_summary=Xsummary;

% contrasts to be calculated:
oat.first_level.contrast={};
oat.first_level.contrast{1}=[3 0 0 0]'; % motorbikes
oat.first_level.contrast{2}=[0 1 1 1]'; % faces
oat.first_level.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
oat.first_level.contrast_name{1}='motorbikes';
oat.first_level.contrast_name{2}='faces';
oat.first_level.contrast_name{3}='faces-motorbikes';
oat.first_level.report.first_level_cons_to_do=[3 1 2];
oat.first_level.post_tf_downsample_factor=4;

oat = osl_check_oat(oat);

%%%%%%%%%%%%%%%%%%%
%% RUN OAT

oat.to_do=[1 1 0 0];
oat = osl_run_oat(oat);

%%%%%%%%%%%%%%%%%%%
%% VIEW OUTPUT

% load first-level GLM result
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

% output freq bins used:
disp('Freq bins in Hz:');
disp(stats.frequencies);

%%%%%%%%%%%%%%%%%%%
%% visualise using Fieldtrip
% note that this produces an interactive figure, with which you can:
% - draw around a set of sensors
% - click in the drawn box to produce a plot of the time series
% - on the time series plot you can draw a time window
% - and click in the window to create a topoplot averaged over that time
% window (which is itself interactive....!)

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.modality='MEGPLANAR';
S2.first_level_contrast=[3];
S2.cfg.colorbar='yes';
S2.cfg.zlim = [-5 5];
S2.view_cope=0;

% calculate t-stat using contrast of absolute value of parameter estimates
[cfg, data]=osl_stats_multiplotTFR(S2);
title([oat.first_level.contrast_name{S2.first_level_contrast}]);

%%%%%%%%%%%%%%%%%%%
%% to do a topoplot averaged over 130 to 160 ms, and 8 to 12 Hz.

cfg.xlim        = [0.13 0.16]; % time window in secs
cfg.ylim        = [8 12]; % freq window in Hz
cfg.interactive = 'no';
figure; ft_topoplotTFR(cfg,data);
title([oat.first_level.contrast_name{S2.first_level_contrast}]);












































%%%%%%%%%%%%%%%%%%%
%% ANSWER/CHEAT: DO SENSOR SPACE SINGLE-BAND TIME-FREQ ANALYSIS USING OAT
% Here we zoom in on single frequency band (5 to 20 Hz) and look at the power of the
% acitivity over time in that band.
% To do this we change S2.num_freqs to be 1, and set the frequency range
% appropriately

oat=[];
oat.source_recon.D_continuous=spm_files;
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
oat.source_recon.D_epoched=spm_files_epoched; % this is passed in so that the bad trials and bad channels can be read out
oat.source_recon.freq_range=[5 20]; % frequency range in Hz
oat.source_recon.time_range=[-0.25 0.5];
oat.source_recon.method='none';
oat.source_recon.dirname=[oat.source_recon.D_epoched{1} '_hilbert_singleband'];

oat.first_level.tf_method='hilbert'; % can be morlet or hilbert
oat.first_level.tf_num_freqs=1; % we are keeping this unusally low in the practical for the sake of speed
oat.first_level.tf_hilbert_freq_res=diff(oat.source_recon.freq_range);
oat.first_level.post_tf_downsample_factor=4;
oat.first_level.time_range=[-0.1 0.3]; % need to make this time range smallet than oat.source_recon.time_range to remove edge effects

% Xsummary is a parsimonious description of the design matrix.
% It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to croat_settingse the (num_regressors x num_trials) design matrix:
Xsummary={};
Xsummary{1}=[1 0 0 0];Xsummary{2}=[0 1 0 0];Xsummary{3}=[0 0 1 0];Xsummary{4}=[0 0 0 1];
oat.first_level.design_matrix_summary=Xsummary;

% contrasts to be calculated:
oat.first_level.contrast={};
oat.first_level.contrast{1}=[3 0 0 0]'; % motorbikes
oat.first_level.contrast{2}=[0 1 1 1]'; % faces
oat.first_level.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
oat.first_level.contrast_name{1}='motorbikes';
oat.first_level.contrast_name{2}='faces';
oat.first_level.contrast_name{3}='faces-motorbikes';

oat.first_level.diagnostic_cons_to_do=[3 1 2];

oat = osl_check_oat(oat);

%% now run the OAT

oat.to_do=[1 1 0 0];
oat = osl_run_oat(oat);

% load GLM result
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

%% visualise using Fieldtrip
% note that this produces an interactive figure, with which you can:
% - draw around a set of sensors
% - click in the drawn box to produce a plot of the time series
% - on the time series plot you can draw a time window
% - and click in the window to create a topoplot averaged over that time
% window (which is itself interactive....!)

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.modality='MEGMAG'; 
S2.first_level_contrast=3;

% calculate t-stat using contrast of absolute value of parameter estimates
[cfg, data]=osl_stats_multiplotER(S2);
title([oat.first_level.contrast_name{S2.first_level_contrast}]);

